C
C     入射平面波在VSH下的展开
C
C  ------------------------------------------------
C  calculating incident wave expansion coefficients
C  for each sphere
c  The p0 and q0 here are the corresponding quantities from
c  xuy95 multiplied by the Mie coefficients an and bn respectively.
c  (modulo normalization), compare xuy95, eqs. 14, 15, 30, 33, 85
c  For x-polarization:
c  p0(n, 1)/an
c  = q0(n, 1)/bn
c  = 1/2 * sqrt(2 n + 1) exp(i k z_l)
c
c  p0(n, -1) = - p0(n, 1)
c  q0(n, -1) = q0(n, 1)
c
c compare eq. 8 of the manual
C  ------------------------------------------------
      subroutine incident_pw(nL, p0, q0, as, bs, aMie, bMie, nmax, k, r0
     $     , indpol)
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2)

      integer nmax(nLp)
      double precision k, r0(6, nLp)
      complex*16 as(nLp, nmp), bs(nLp, nmp), p0(nLp, nmp), q0(nLp, nmp),
     $     aMie(nLp,np), bMie(nLp,np),A,ci,cin,cmz

      common/fnr/    fnr(0:2*(np+2))

      ci     = dcmplx(0.d0, 1.d0)
      cin    = dcmplx(0.d0,-1.d0)
 18   do imn = 1,nmp
         do i = 1,nL
            p0(i,imn) = 0.d0
            q0(i,imn) = 0.d0
            as(i,imn) = 0.d0
            bs(i,imn) = 0.d0
         enddo
      enddo

      do 20 i=1,nL
         cz = dcos(k*r0(3,i))
         sz = dsin(k*r0(3,i))
         cmz= 0.5d0*dcmplx(cz,sz)
         do 21 n=1,nmax(i)

C     the index of a_{l,1} and b_{l,1}
            imn = n*n + n + 1

            A = fnr(2*n+1)*cmz
C
C     for each order n, only the -1 and 1 index are nonzero
C     and it is the same for all spheres
C
            p0(i,imn) = aMie(i,n)*A
            q0(i,imn) = bMie(i,n)*A
            p0(i,imn-2) = -p0(i,imn)
            q0(i,imn-2) =  q0(i,imn)

            if(indpol.gt.1) then
               p0(i,imn) = p0(i,imn)*cin
               q0(i,imn) = q0(i,imn)*cin
               p0(i,imn-2) =  p0(i,imn)
               q0(i,imn-2) = -q0(i,imn)
            endif

C     copy the incident wave as initial guess
            as(i,imn) = p0(i,imn)
            bs(i,imn) = q0(i,imn)
            as(i,imn-2) = p0(i,imn-2)
            bs(i,imn-2) = q0(i,imn-2)
 21      continue
 20   continue


      return
      end
